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ABSTRACT 


A description of the FNWC five level baroclinic, global, primitive 
equation model is presented including the finite difference equations 
used. Analytic initial data is generated in order to avoid the compli- 
cations caused by real data with respect to verification. A comparison 
of second order differencing versus mixed second and fourth order 
differencing shows that the latter gives relatively more accurate 
propagation for smaller scale disturbances than the former. In addi- 
tion various experiments are conducted to improve the model's treatment 


of cross-polar flow. 
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I. INTRODUCTION 


After the FNWC Northern Hemispheric Primitive-Equation Model became 
operational (Kesel and Winninghoff, 1972), Dr. F. J. Winninghoff, a 
former Naval Postgraduate School faculty member, designed and programmed 
a five-level global primitive equation model utilizing a spherical, 
staggered grid and the sigma coordinate system. This global atmospheric 
model was basically a conversion of the operational hemispheric model 
to a global model using a staggered latitude-longitude grid instead of 
a polar stereographic projection. 

The advantage of a global model is that it treats the entire atmos- 
phere as a physical system. It is reasonable to expect that when the 
model is fully operational, increased accuracy and range of prediction 
will result. A global atmospheric model allows tropical circulation 
to develop freely without artificial boundary conditions. This is crucial 
since the tropics are the main source region of the energy, which drives 
the general circulation of the atmosphere. Interactions between the 
hemispheres are also represented. 

The operational use of a global atmospheric model for short range 
prediction is practical since advanced computers capable of integrating 
a global model in a reasonable amount of time are presently available 
and since satellites in the near future are expected to provide more 
accurate world wide data input. 

From a tactical standpoint an operational global atmospheric model 


is a tremendous advantage since it provides a timely forecast of the 
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weather anywhere in the world. Improved predictability in the tropics 
would greatly benefit the Navy since a large percentage of its bases 

are located in the tropics and many naval operations take place in 
tropical waters. The most recent example of this is the establishment 
of U. S. Naval presence in the Indian Ocean and the building of a naval 
base on Diego Garcia Island in that area. 

Elias (1973) modified and carried out further experiments with the 
Winninghoff global model on the 6500CDC computer located at FNWC. He 
generated geopotential fields by inserting an analytic spherical harmonic 
stream function (Neamtan, 1946) into the linear balance equation and 
solving for the geopotential by successive over-relaxation. The geo- 
potential fields were generated on a 63 x 63 polar stereographic grid 
and were interpolated onto a 5° Northern Hemisphere latitude—longitude 
grid after insertion into the global model. Once this interpolation 
had taken place the fields were reflected into the Southern Hemisphere. 
Real data analyzed by FNWC objective schemes were also inserted into 
the global model on a polar stereographic grid and then interpolated 
onto a spherical grid and reflected into the Southern Hemisphere. In 
the cases where analytic data were generated,artificial moisture and 
temperature fields were also produced. The purpose of using an analytic 
geopotential field instead of real data is that it allows pa number, 
phase speed and wave amplitude to be specified. 

The geopotential fields were used in the linear balance equation to 
obtain the stream function which provided initial rotational winds. 
In the analytic cases Elias also avoided numerical balancing by using 


an analytic solution for the wind fields. His forecasts using analytically 
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derived winds were well behaved but those using winds from the linear 
balance equation excited spurious inertial-gravity waves which were 
undesirable for operational forecasts, 

Mihok (1974) and McCollough (1974) converted the FNWC global atmos- 
pheric model to operate on the IBM 360 computer at the Naval Postgradu- 
ate School. McCollough examined various schemes of initialization of 
the model using real data from the FNWC objective analyses used by 
Elias. He found that the use of the Robert (1965) time frequency fil- 
ter and dynamic balancing using forward and backward averages about 
the initial value was effective in reducing inertial-gravity wave 
"noise." 

McCollough (1974) pointed out that the method of data input by 
interpolating to latitude-longitude grid points and to sigma surfaces 
from a polar stereographic grid on p-surfaces probably introduces some 
imbalance between mass and wind fields. In order to avoid the 
problems associated with initializing with real data this study employs 
a simplified direct analytic solution to the non-linear balance 
equation (Phillips, 1959). This solution is based on the same spheri- 
cal harmonic stream function (Haurwitz, 1940, Neamtan, 1946) used by 
Elias (1973). 

McCollough also noted that large gradients developed near the 
pole when real data were used. This thesis examines the polar problem 
with an analytic initial state which gives flow over the pole. 

Mihok (1974) modified the horizontal spatial difference equations 
to form a mixed second and fourth order scheme with the aim of improving 


the phase speeds of the meteorological waves, as proposed by Williams 
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(1972). He generated analytic fields using the same method as Elias 
(1974) and also used the same real data as Elias (1973) and McCollough 
(1974). Mihok derived his winds from the stream function computed by 
numerically solving the nonlinear balance equation with geostrophic 
winds used in the Jacobian. The use of this method of balancing intro- 
duced inertial-gravity waves which made it difficult to compare the 
relative merits of forecasts made with the second and fourth order 
scheme. The experiments in this thesis avoid the balancing problem by 
employing analytic height and wind fields. Mihok's mixed scheme is 


compared with the second order scheme for various wave numbers. 
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II. BAROCLINIC PRIMITIVE EQUATION MODEL 


The five-level global atmospheric primitive equation model currently 
under development by FNWC and the Naval Postgraduate School is used in 
this study. The primitive equations are written in spherical coordinates 
with sigma used vertically. They are similar to the equations used by 
Siiiorinsiy et al (1965), Arakawa et al (1969) and Kesel and Winninghoff 
(1972). Finite differencing of these equations takes place on a grid 
staggered both horizontally and vertically. 

The finite difference form of the primitive equations contain Ara- 
kawa's (1966) technique for conserving energy. Fictitious kinetic 
energy is not produced by the integration of the nonlinear advective 
terms. However, difference form of the hydrostatic equation used in 
the model does not meet the restrictions placed on differencing in the 
vertical if total energy is to be conserved. A discussion of the re- 
quirements for total energy conservation may be found in Haltiner (1971). 
The difference equation used in the model is similar to the one used 
by Kesel and Winninghoff (1972). They state that the reason for 
differencing the hydrostatic equation in this manner is to reconcile 
the initial temperature, He ct and surface pressure analyses with the 
forecast variables in the model. 

A complete list of the finite difference equations, similar to the 
list given by Mihok (1974), is provided in Appendix A. The fourth 
order modifications to the difference equations (Mihok, 1974) are stated 


in Appendix B. 
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Tke heating and moisture source terms, lateral diffusion terms and 
frictional stress terms are not included in Appendix A but are the 
same as those used in the operational FNWC Northern Hemispheric Primi- 
tive-Equation Model (Kesel and Winninghoff, 1972). For the purpose of 
this study the heating and moisture source terms and lateral diffusion 
terms are turned off. 

A flat earth (all terrain heights at sea level) is used in this 
study since it was desirable to avoid the instability experienced by 


McCollough (1974) when terrain height was included. 


A. PRIMITIVE EQUATIONS 


The continuous form of the equations used in the model are as follows: 


Longitudinal momentum equation 


(ut) _ _ 1 jue + d(uvT cos @)) 
ot a cos 6 oA ele) 


To (wu) 4 Muv_tan 8 


+ AG + Tvft 
1 109 oT 
“a cos 0 [ oA pee an! = ue ie Ye (1) 
Latitudinal momentum equation 
(vm) _ _ if pocuvT) ¥ d(vvT cos 9); 
ot a cos @ oA 00 
et Td(wv) _ Tuu tan o teak 
00 a 
a7 36) on 
aie [1 6 + RT 76! + TF + dD, (2) 
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Thermodynamic equation 


ont 1 d(muT) , I(mvT cos 6) 


= [ 1d (wT) 
ot a cos 8 oA 00 


as 00 


RT oT 1 oT oT 
trang [tet Ge Sais tea" [u 5) + v 39 cos 91)] 
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Moisture continuity equation 





oqT 5 AGED + o(tvq); 
dt a cos 30 


+ e + 7Q + Dq 


Mass continuity equation 


on” 1 0 (uT) 


- + 


an [ d(vT cos 8) ow 
ot a cos 8 oA 00 


J+ AG 


Hydrostatic equation 


cl’) a BE 


30 o 
Equation of state 
aP = RT 
The dimensionless vertical coordinate sigma (co) is defined as 


Sey ta 
T 


The measure of the vertical motion 
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(4) 


(5) 


(6) 


(7) 


(8) 
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A complete list of symbols and abbreviations used in the above equations 
may be found in the front of this thesis. 

The terrain pressure derivatives appearing as part of the pressure 
gradient force terms in the momentum equations are computed by a special 
finite difference scheme in order to reduce the truncation error noted 
by Kurihara (1968). Even though the main purpose of this scheme is to 
correct problems occurring over steep terrain, it was allowed to remain 
in the model during the experiments reported in this thesis in which a 
flat earth (all terrain heights at sea level) was specified. The tech- 
nique involves the expression of the gradient of terrain pressure as a 
function of geopotential on pressure and sigma surfaces. Equation (10) 


illustrates this relationship. 


1 1 om om T 
= — f. —i —j = — — 
Vn a ‘cos 8 dA~ ‘si 502! RT [Vo VO, ete? 





The subscript p represents a pressure surface and the subscript oO 
represents a oO surface. The finite difference scheme given in Appendix 
A utilizes this relationship by interpolating the value of the geo- 
potential at adjacent grid points to the pressure surface of the com- 
putational grid point and then performing the integration on a pressure 
surface. The computational grid point is common to a sigma surface and 
the pressure surface on which the integration takes place. Thus, pres- 
sure force is computed at a point on the sigma surface, 

The terrain pressure derivatives appearing in the thermodynamic 
equation as part of the term representing the work done per unit time 
and mass by the pressure force are also handled in a special way as 


shown in Appendix A. 
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B. SPHERICAL COORDINATES, STAGGERED GRID 

The spherical grid points occur at five degree intervals of latitude 
and longitude. These grid points are subdivided into two categories; 
mass points at which the mass variables T, 9, T, q, dad w are retained 
and velocity points at which the motion variables u and v are carried. 
The poles are considered mass points. Mass points and velocity points 
are staggered throughout the grid, hence mass and velocity points occur 
at ten-degree intervals. Figure (1) taken from Mihok (1974) illustrates 
the spherical coordinate, staggered grid. 

When a variable is needed at a grid point on which it is not 
carried, an average of the values at the surrounding grid points is 
taken. For example if 1 is needed at a velocity point, the average 
of the tm values at the four mass points surrounding the velocity point 


is used. 


C. VERTICAL LAYERS 
The atmosphere is vertically divided into five layers in sigma. 

The variables u, v, T, > are retained at the 0.9, 0.7, 0.5, 0.3 and 0.1 
sigma levels, which are the centers of sled of the five layers. The 
specific humidity (q) is retained only in the lower three layers. This 
is reasonable for the distribution of moisture in the atmosphere. The 
quantity w=- G6 is carried at the interfaces of the vertical layers 
(0.8, 0.6, 0.4, 0.2 sigma levels) and is zero at the upper and lower 
boundaries. Figure 2 provides a graphic description of the vertical 


atmospheric layering used in the model. 
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D. SPATIAL FINITE DIFFERENCE METHODS 

The experiments in this thesis fall into two categories: (1) experi- 
ments comparing second order spatial differencing versus mixed second 
and fourth order spatial differencing, and (2) experiments examining the 
reaction of the model to flow over the poles with second order differenc- 
ing. 

1. Second Order and Mixed Second and Fourth Order Finite Differencing 

Williams (1972) has shown that phase speed errors with fourth 
order finite differencing are considerably less than phase speed errors 
occurring with second order differencing. However, fourth order differ- 
encing requires a reduction in the length of the time step in order to 
insure linear computational stability. (Time step requirements are dis- 
cussed more fully in Section II-E.) Williams found that the time step 
would not have to be reduced as much if fourth order differencing is 
limited to the advection and divergence terms. These terms are the prin- 
cipal ones which affect the phase speeds of the meteorological waves. 

In-order to gain improved phase speeds for meteorological waves 
while keeping the necessary reduction of the time step to a minimun, 
Mihok (1974) inserted a mixed second and fourth order differencing scheme 
into the model. The mixed scheme is documented in Appendix B. 

Another factor affecting the phase speeds of waves in the model 
is the. arrangement of points in the horizontal. Phase speed errors for 
waves moving along the coordinate axes, with the difference scheme of 
this thesis, are the same as those for an unstaggered grid. However, 


Williams (1972) has noted that waves propagating at an angle of 45° to 
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the coordinate axis will have a larger phase speed error on the staggered 
grid than on the unstaggered grid. This is because the staggered grid 
points used in the difference equations are more widely separated rela- 
tive to a wave that is moving at a 45° angle with respect to the coordi- 
mate axes. 

2. Polar Finite Differencing 

One of the difficult problems associated with spherical coordinates 
is the treatment of the poles which are singular points. The zonal and 
meridional wind components are undefined at the poles. In the difference 
equations for this model any term requiring a velocity component from the 
poles is set to zero. 

Since the poles are considered mass points in this model separate 
calculations of the mass point parameters (1, T, q, $, w) are made for 
the poles. Following the procedure of Arakawa (1974), polar mass point 
parameters in the model change as the result of the meridional mass flux 
at all velocity points on the 85° latitude circle. The flux at each mass 
point is weighted to represent a triangular area with apexes located at 
the pole and at the mass points adjoining the velocity point. When an 
individual mass point parameter has been calculated for each of the 36 
triangular areas surrounding the pole, an average is taken and inserted 
as the pole value. As an illustration, the equation for the contribution 
of one atmospheric layer to the local change in terrain pressure can be 


written 


1 Aut , A(vt cos 0) 


a os ae Aut 
ot a eae SAK SAE 4 (11) 


where <> indicates a vertical summation of the layer contributions 
and 7 the terrain pressure at a velocity point obtained by averaging 


the 1 values at the four surrounding mass points. 
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By clearing all terms from the denominator of (11) and multiply- 


ing by the earth's radius, a , one can obtain 


(32) [a” 4AAA@ cos 6] = < Aut a2h0 + A(vii cos 8)a2Md > — (12) 


The quantity on the left hand side, ee 4AXAA8 cos 8], represents the 
approximate area of a ten-degree square with a velocity point at its 
center. The u and v values at the central velocity point control the 
flux over the area. The triangular flux areas surrounding the pole are 
one quarter the size of the square flux areas at lower latitudes. Thus, 


the equation for each polar triangular flux area becomes 


(2) [a? AAAe cos 8] = [+< vil cos > a” AA] oa 


All mass flux terms except the meridional flux along the 85° latitude 
circle are omitted . 

The following equation shows how the contribution of each of 
the 36 flux triangles is averaged to obtain the final value for the 


local change in terrain pressure at the pole. 


36 

ene IeaQuse DL! = : 

3 aon a ge ae A (14) 
n=1 


An experiment was also performed using values at the poles ob- 
tained solely from averaging the desired quantity around the 85° latitude 
Keekon This method {2 much simpler to apply and requires less computer 
time. 

Another problem which occurs near the polar region is maintaining 


computational stability as the meridians converge to the poles and the 
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grid distance along a latitude circle steadily decreases. Since it is 
not desirable to reduce the time step or modify the grid Arakawa's (1974) 
method of smoothing the zonal component of the pressure gradient force 
and the divergence is used. A stability coefficient, given in the 


following Equation (15) is used to determine the smoothing. 


D cos 8 
Bn C, At sin (md) (15) 
where 
S = stability coefficient 
D = 5° grid distance at equator ) 


C_ = phase speed of the fastest gravity wave 
At = time step 
m = wave number 


d = 5° in radians 


If the value of S is less than 1 smoothing is necessary. The field 
to be smoothed is expanded into a Fourier series and the amplitude of 
each unstable wave component is reduced by a factor proportional to S. 
Arakawa states that this method of smoothing does not smooth the fields 
of variables because it is simply a generator of multiple point differ- 
ence quotients. 

Several of the experiments with flow over the pole involved 


modification of the stability coefficient or direct smoothing of the main 


variables in addition to the zonal pressure gradient and divergence fields. 
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E. TIME FINITE DIFFERENCE METHODS 
The initial time step in the model is performed using the Euler back- 


ward (Matsuno) difference scheme which is expressible as follows: 


* t gre 
F =F +At Set [first approximation] 
; e é (16) 
ai = F° + At [forecast] . 


F is a vector including the dependent variable, the superscript t 
denotes the time step, the superscript * indicates the result of the 
intermediate step and At is the time interval of a single time step. 

The Euler backward scheme is a two step iterative scheme which is either 
damping or neutral depending on the wave number (Haltiner and Williams, 
1973). The principal damping takes place in short waves. For this 
reason the Euler backward scheme is applied at three-hour intervals in 
addition to the initial application. 

All other time steps are performed using central differences, commonly 


called the leapfrog scheme; 


itl 2 pols ont a Sale iak) 
The leapfrog scheme gives rise to a spurious wave called the computational 
mode which has no counterpart in the analytic solution to the differ- 
ential equation. The phase of the computational mode changes every time 
step causing the solutions at even and odd time steps to decouple. The 


amplitude of the computational mode is greatest for short wavelengths 


(Haltiner, 1971). 
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In order to limit the effects of the computational mode and reduce 
noise in the short wavelengths the Euler backward scheme is applied 
every three hours. As previously mentioned, this tends to dampen short 
waves somewhat. In addition, it is important to note that there is 
no computational mode associated with the Euler backward scheme and 
therefore application of this scheme recouples the solution. 

A time step of ten minutes is used when second order differencing 
is applied throughout the model. When the mixed second and fourth 
order differencing scheme (Williams, 1972), (Mihok, 1974) is used the 
time step is reduced to eight minutes. Williams showed, experimenting 
with the linearized barotropic primitive equations, that a 14% smaller 
time step was needed with mixed second and fourth order spatial 
differencing in order to meet the von Neumann necessary condition for 
stability. Both the ten-minute and eight-minute time steps are too 
large near the poles to maintain stability unless Arakawa's method of 


smoothing derivatives is used. 
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III. INITIALIZATION 


Initial conditions for all experiments in this thesis were derived 
analytically. They approximate an initial barotropic state. Generation 
of initial conditions in this study differs from the method used by 
Elias (1973) and Mihok (1974) in that an analytic solution for @ rather 
than a numerical solution of the balance equation is used. Generation 
of initial conditions in this manner was found to be less time consuming 
and produced waves with almost no noise that propagated with small dis- 
tortion. 

Initial temperatures are generated on twelve constant pressure 
levels distributed from 1000 mb to 50 mb and height values are generated 
for ten of these levels. In addition, sea level pressure and sea surface 
temperature fields are produced. Wind fields are generated at seven 
constant pressure levels. All generated fields are interpolated to 
sigma surfaces. Initial values at the poles are averages taken around 


the 85° latitude circle on sigma surfaces. 


A. ANALYTIC BALANCING 

The initial geopotential and velocity fields used for the experiments 
in this thesis were obtained from the stream function solution to the 
linearized vorticity equation (Haurwitz, 1940).Neamtan (1946) later 
solved the complete vorticity equation and obtained results verifying 


Haurwitz's work. 
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The expression for the stream function, ~ , is given by 
* 2 
Wy =A _ sin (m-vt) sin 0 cos m6 - B a sin 8 (18) 


where m is the wave number, vy is the angular wave velocity, a is 
* 
the radius of the earth and A and B are constants. The amplitude 


* 
constant, A , (Elias, 1973) is defined as 
* N 
A =A [(2m!/2° N!)] , (19) 


where N=m+l and A is aconstant. The constant B is directly 


related to the angular phase speed by the equation 


Ve p NG) - 2 22 


n N(N+1) ~ N(N#1) (20) 


where (v/m) is the angular phase speed, and 2 is the angular velocity 
of the earth. 

Haurwitz (1940) has shown that harmonic waves defined by this stream 
function in a barotropic non-divergent atmosphere will move with a 
constant angular velocity without changing shape. Harmonic waves 
generated in this way Srovide excellent controlled conditions for testing 
and debugging atmospheric models. Many experimenters including Phillips 
(1959), Gates (1962), Heburn (1972), Holloway, Spelman,and Manabe (1973), 
Elias (1973), and Mihok (1974) have used the Haurwitz stream function to 
generate initial conditions. 

It should be noted that the equations used in the model are not those 
of a non-divergent atmosphere. Phillips (1959) has noted that the 


presence of divergence in the barotropic atmosphere will slow the rate 
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of progression of the flow pattern especially for small values of the 
wave number m. 

The Haurwitz stream function is used by Phillips (1959) as the 
forcing function to obtain the geopotential from the non-linear balance 
equation 


V4 = ey - We ve - 20)? 4 2A BY (21) 
dx? dy? 


Equation (21) in spherical coordinates becomes 











1 1 926 P) od 
5S =q (cos'6)-~)] = 
az cos26 92 cos 6 300 00 
f£ 1 o7y Levu ow 1 op Of 
oF eS se (cos 8) | FP oe ee 
Bec tp (gy? , cos 6 ae 30 a2 00 9 
1 o7y .2 2 o7y a2p 
- 2065 Ea ee es ee nee? 
a 6 a cos @ 


The solution for 9 of the non-linear balance equation may be stated as 


ie a K(0) + a HCO) sin m+ 2 -C(O) (2 sin” mi = 1) (23) 


where @$ is the geopotential perturbation 


2 
ot GCG )con 8 + (n-=m—2)= as ay 
cos 0 


i re 2 
+4 





A(®) = 5 (20+ Bycos- 


(24) 
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2(0+B) AL 
a" 4 2 mb. phiny 
B(6) = (atl) Gnt2) cos 9 [(m + 2m+#2) - (mtl1)~ cos’ 6] (25) 
and 
ew 2 
c(é) =< 4? cae ae biliGel coat ke) = vata (26) 
a 


This solution was obtained by Phillips (1959). Equation (23) provides 
a geopotential perturbation field. The height at a particular pressure 
level is obtained by adding the height perturbation from (23) to the 
mean height for that level obtained from the NACA standard atmosphere 
(Haltiner and Martin, 1957). 

The linear balance equation may be obtained by dropping the last 
two terms in (21) for Cartesian coordinates and (22) for spherical co- 


ordinates. The linear solution is given by 


@ = F(6) + G(8) sin (m - vt) (27) 
where 
2 2 
F(6) = BRa cos 8 (28) 
and 
& m* + 2m + 2 m 
G(8)sin(mA-vt) =A 221) (m2) sin(mA-vt) cos 9 
Ps ao Ao sin (m\ - vt) cast ) (29) 
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The solution to the linear balance equation produces initial fields 
which are in geostrophic equilibrium while the non-linear solution pro- 
duces fields which are in approximate gradient equilibrium and therefore 
more realistic. This is especially true near the equator where the 
geostrophic approximation is invalid. For this reason non-linear 
balancing is used in all experiments in this thesis except Experiment 3 
which involves a flow pattern of wave number 12. 

In Experiment 3,non-linear balancing results in aliasing, (Haltiner, 
1971), which distorts the initial fields. This problem was eliminated 
by the use of linear balancing. The effect of eliminating the non- 
linear terms is minimal because the wave amplitude used in the experi- 


ments is small. 


B. ANALYTIC WINDS 
Analytic zonal and meridional wind components are obtained from the 


derivatives 


sat 
gia ae (30) 


and 


1 ow 


acos 0 oA : (31) 


The results are 


* ms 
- 2[a* sin(m-vt) cos™' "6 —mA sin(mA-vt)cos” 15 51n26-Ba~cos0] 


c 
" 


(32) 


< 
" 


2fa"n sind eee 8 cos(mA - vt)] . (33) 
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C. ANALYTIC TEMPERATURES 
Analytic temperatures, constant on each pressure surface, are 
generated in accordance with the National Advisory Committee for 


Aeronautics (NACA) standard atmosphere (Haltiner and Martin, 1957). 
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IV. RESULTS 


The experiments conducted for this thesis fall into two categories: 
(1) experiments comparing second order differencing versus mixed 
second and fourth order differencing as discussed in Section LI-D-1 
and (2) experiments attempting to improve the forecast effectiveness 
of the model when there is flow over the poles. 

In all experiments presented in this thesis analytically derived 
fields are used as described in Section III. The variation in the 
character of the initial harmonic waves was accomplished by varying 
the specification of the wave number, phase speed and the amplitude 
constant A which is part of the coefficient Ti given in Equation (19). 

The use of analytic fields resulted in the elimination of Ps 
inertial-gravity waves resulting from initial inconsistencies between 
the wW and 9 fields. Thus, results directly related to the finite 
differencing of the primitive equations can be examined under con- 
trolled conditions without distortion due to initialization. 

Surface pressure analysis and forecasts are used to illustrate the 
output of the model. These analyses,for the Northern Hemisphere only, 
were obtained by interpolating the spherical grid onto the standard 
63 x 63 polar stereographic grid used by FNWC. The interpolation of 
the spherical data field for a point on the 63 x 63 grid is accomplished 
by using the Ayres’ central difference formula, which is continuous in 
the first derivative. If the point lies within the border zone, a 
double linear interpolation is performed. Very little detail is lost, 


except at the pole as shall be shown in Experiment VIII. 
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Harmonic analysis (Jeffreys and Jeffreys, 1956) is performed on the 
surface pressure analysis and on the forecasts at three-hour intervals. 
The harmonic analysis reveals the wave numbers present around a latitude 
circle and their amplitudes and phases. 

A. SECOND ORDER DIFFERENCING VERSUS MIXED SECOND AND 

FOURTH ORDER DIFFERENCING 

Experiment I. In this experiment,wave number 6, a phase speed of 
+ 12° longitude per day and A = 1000 are used. Chart A shows the initial 
surface pressure field. Chart B is the 48-hour forecast using second 
order differencing and Chart C shows the 48-hour forecast using the mixed 
differencing scheme. Note that the forecast fields retain their defini- 
tion and are not distorted by undesirable inertial gravity waves. 

Figure 3 shows the phase angles in degrees longitude as a function of 
latitude at 12-hour intervals to 48 hours for second order differencing. 
Figure 4 shows the same relationship for the mixed differencing scheme. 
Comparison of Figures 3 and 4 shows a considerable improvement in the 
forecast phase speed of the waves when the mixed differencing scheme is 
used. This is especially true as the latitude increases. Above 60° 
latitude the mixed scheme forecasts of phase speed are too fast while 
the second order forecasts indicate a greatly retarded phase speed. 

Experiment II. Wave number 9, a phase speed of + 12° longitude per 
day and A = 0.1 are used in this experiment. Chart D shows the initial 
surface pressure field. Chart E is the 48-hour forecast using second 
order differencing. Chart F is the 48-hour forecast using the mixed 


differencing scheme. The forecast and initial fields are not distorted 
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by undesirable gravity waves. Figure 5 shows the phase angles in degrees 
longitude as a function of latitude at 12-hour intervals out to 48 hours 
for second order differencing. Figure 6 shows the same relationship for 
the mixed differencing scheme. Comparison of Figures 5 and 6 indicates 

a closer approximation to the true phase speed when the mixed scheme is 
used. Above 60° latitude the mixed scheme phase speeds are too fast while 
considerable slowing is evident for second order differencing. 

Experiment III. Wave number 12, a phase speed of + 12° longitude per 
day and A=5x 1000 are used in this experiment. The initial surface 
pressure field is shown in Chart G while Charts H and I are the 48-hour 
forecast fields for the second order differencing and mixed differencing 
schemes, respectively. Note that the initial and forecast fields are 
not distorted by gravity waves. Figure 7 shows the phase angles in 
degrees longitude as a function of latitude at 12-hour intervals out to 
48 hours for second order differencing while Figure 8 shows the same 
relationship for the mixed differencing scheme. When the second order 
differencing is used there is very little wave propagation. The phase 
progression is very distorted and marked retrogression is noted above 
40° latitude. Figure 8 indicates a distinct improvement in the propa- 
gation of the waves when the mixed differencing scheme is employed. In 
Figure 8 the phase speed of the waves increases with latitude. However, 
at all latitudes the phase speeds are slower than for the nondivergent 
analytic solution. 

It should be noted that the solution to the linear balance equation 
was used to generate the initial height fields in this case. The use 


of the non-linear solution resulted in undesirable aliasing. 
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B. CROSS-POLAR FLOW 

McCollough (1974) remarked that "noise" developed in the polar 
regions as the forecast time progressed and that longer integrations 
indicated a need for a review of the polar region finite differencing 
and the method of handling the pole point. In order to examine the 
effectiveness of the model in forecasting for the polar regions, flow 
over the poles was generated analytically using wave number 1 and a 
phase speed of -94.4° longitude per day. After some experimentation 
it was decided to make A = 3 x 107 yielding a maximum wind of 20 kts 
in the polar region. This method of generating cross-—polar flow was 
taken from Holloway, Spelman and Manabe (1973). 

Chart J is the initial surface pressure analysis used for the 
flow over the pole experiments. Charts K-N illustrate surface pressure 
forecasts at 12-hour intervals out to 48 hours. All cross polar flow 
experiments use second order differencing. 

Examination of Chart J reveals that the isobars are evenly spaced 
and curve smoothly over the pole. However, the pattern of the isobars 
becomes increasingly distorted near the pole with progressive fore- 
casts as can be seen in Charts K-N. The objective of the subsequent 
experiments is to test various modifications to the model in order to 
improve the accuracy of carcerare in the polar region when cross-polar 
flow is present. 

Experiment IV. Based on the assumption that the distorted flow 
pattern in the polar regions was the result of errors in the derivatives 
near the pole, smoothing of the derivatives was increased by squaring 


the stability coefficient, S , defined in Equation (15). This 
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modification did not result in smoother forecasts in the polar regions 
and in fact the disturbance was slightly augmented as can be seen by a 
comparison of Chart L and Chart 0. 

Experiment V. In this experiment removal of the irregularities in 
the polar region was attempted by applying the Arakawa smoothing tech- 
nique directly to the main variables at three-hour intervals. The 
Arakawa smoothing was applied to T, u, v, T and q on the ninth and 
tenth time steps after the Matsuno step. Charts P and Q are the 24-hour 
and 48-hour forecasts made with this technique. Comparison of the above 
charts with Charts L and N reveal a reduction in the amplitude of the 
disturbance at the pole and a more uniform flow at the higher latitudes. 

Experiment VI. The promising results of Experiment V made it seem 
probable that application of Arakawa smoothing to the main variables 
at every time step would totally remove the distortion of the isobars 
near the pole. Although this method gave somewhat improved results over 
the case where the main variables were not smoothed, the effect was 
still not as good as that achieved in Experiment V. 

Experiment VII. Holloway, Spelman and Manabe (1973) transformed 
vector components into a polar stereographic mapping prior to Fourier 
filtering in order to prevent excessive smoothing of flow over the pole. 
Their Fourier filtering method is similar to the Arakawa smoothing tech- 
nique except that the shorter waves are completely removed. In this 
experiment the main variables were smoothed every three hours, as in 
Experiment V; but prior to smoothing, the wind components u and v were 
transformed to a polar stereographic projection plane. The method of 


transformation is the same as used by Holloway, Spelman and Manabe. 
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This experiment resulted in the model becoming unstable after 27 
hours. 

Experiment VIII. At this point in the experimentation it was noted 
that the interpolation technique used to transform the surface pressure 
field from spherical coordinates to a 63 x 63 polar stereographic grid 
was not adequately handling the pole point. Hence, the interpolation 
scheme was bypassed at the pole and the pole value on the spherical 
grid was inserted directly onto the 63 x 63 grid which is centered at 
the pole. Chart R and S show the result of this modification of the 
output interpolation scheme. The flow appears to be smoother indicating 
that a portion of the distortion in the polar regions was due to a 
faulty output interpolation scheme. 

Experiment IX. This experiment tests the effect of using only 
averages around the 85° latitude circle for values needed at the poles. 
The modified output interpolation scheme (Experiment VIII) is included 
in the model for this experiment. The resulting 24-hour forecast is 
shown in Chart T which shows a considerable improvement over Chart L 
(the original cross-polar flow 24-hour forecast). However, there is 
little difference between the 24-hour forecast using just the modified 
output interpolation scheme (Experiment VIII, Chart S) and Chart T which 
includes the additional feature of using averages for polar parameters, 

Experiment X. This experiment consisted of a combination of the 
modified output interpolation scheme (Experiment VIII) and the applica- 
tion of Arakawa smoothing to the main variables at three-hour intervals 
(Experiment V). Chart U shows that Arakawa smoothing of the main 
variables at three-hour-intervals reduces the distortion at the pole by 


a considerable amount. The improvement was larger than in Experiment V. 
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A 48-hour forecast was not produced in this experiment due to time 


limitations. 
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3 is the column counter and 4 is the row counter. 


Peunt (i,1) is at 85S, 10E. 
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FIGURE 2. VERTICAL LAYERING 


In the above figure, sigma (0) is the dimensionless vertical 
cocrdinate, u and v are the zonal and meridional wind 
components, respectively, q is the specific humidity, $¢ is 
the geopotential, 1 is the terrain pressure and w is the 
vertical velocity (-c). 
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Figure 3. Phase angle (degrees longitude) vs latitude for second order 
differencing, wave number 6, phase speed 12°/day, A = 1000. (Latitudes 
with zero wave amplitude are not included.) 
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Figure 4. Phase angle (degrees longitude) vs latitude for mixed second 
and fourth order differencing for wave number 6, phase speed 12°/day, 
A = 1000. (Latitudes with zero wave amplitude are not included.) 
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Figure 5. Phase angle (degrees longitude) vs latitude for second order 
differencing, wave number 9, phase speed 12°/day, A = 0.1. (Latitudes 
with zero wave amplitude are not included.) 
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Figure 6. Phase angle (degrees longitude) vs latitude for mixed second 
and fourth order differencing, wave number 9, phase speed 12°/day, 
A= 0.1. (Latitudes with zero wave amplitude are not included.) 
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Figure 7. Phase angle (degrees longitude) vs latitude for second order 
differencing, wave number 12, phase speed 12°/day, A= 5 x Oma % 
(Latitude with zero wave amplitudes are not included.) 
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Figure 8. Phase angle (degrees longitude) vs latitude for mixed second 
and fourth order differencing, wave number 12, phase speed 12°/day, 
A=5x 10% . (Latitudes with zero wave amplitudes are not included.) 
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Chart C. 48-hour surface pressure forecast, mixed second and fourth 
order differencing, wave number 6, phase speed 12°/day, A = 1000. 
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‘Chart H. 48-hour surface pressure forecast, second order differencing, 
wave number 12, phase speed 12°/day, A = 5.0 x 2075.2 


56 

















oc. 


ne 


ae 


Pa 


er 






































Chart 0. 24-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A - 3.0 x 10’; Arakawa smoothing coefficient squared. 
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Chart P. 24-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A = 3.0 x 107; Arakawa smoothing applied to main 
variables at 3-hour intervals. : . 
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Chart Q. 48-hour surface pressure forecast, wave number 1, phase 
speed -94,4°/day, A = 3.0 x 10’, Arakawa smoothing applied to main 
variables at 3-hour intervals. 
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‘Chart R. 12-hour surface pressure forecast, wave 
speed -94.4°/day, A = 3.0 x 10’, modified output interpolation 
scheme, 
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24-hour surface pressure forecast, wave number 1, phase 


speed -94,4°/day, A = 3.0 x 10’, 


Chart S. 


modified output interpolation scheme. 
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‘Chart U. 24-hour surface pressure forecast, wave number 1, phase 
speed -94,4°/day, A = 3.0 x 10°, modified output interpolation 


scheme, Arakawa smoothing applied to main variables at 3-hour 
intervals. 
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V. CONCLUSIONS 


The experiments comparing second order differencing and the mixed 
second and fourth order differencing scheme indicate that in general 
the mixed differencing scheme provides a much closer approximation to 
the true phase speed. Phase speeds using second order differencing 
tended to decrease and therefore become more inaccurate with increas- 
ing latitude. The mixed differencing scheme produced phase speeds 
which tended to increase with latitude therefore becoming more accurate. 
As the wavelength of the harmonic waves became shorter the mixed 
differencing scheme gave increasingly better phase speeds in com- 
parison with the second order phase speeds. This is especially evident 
in Experiment III where wave number 12 was used. It should be noted 
however that both schemes when compared only with themselves on the 
basis of wave number became steadily less accurate as the wave number 
increased. 

The distortion in the cross-polar flow was in part due to the 
interpolation scheme which transformed the spherical grid to a 63 x 63 
polar stereographic grid. A great deal of the remaining "noise" 
at the pole can be removed by applying Arakawa smoothing at three- 
hour intervals to the main variables. 

The application of Arakawa smoothing on a polar stereographic plane 
projection for the main vector variables caused the model to become 
unstable. Squaring the Arakawa smoothing coefficient did not have a 


beneficial effect on the cross-polar flow. Application of the Arakawa 
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smoothing technique to the main variables at every time step smoothed the 


cross-polar flow somewhat but was not as effective as the three-hour 


application. The use of 85° latitude averages for the polar quantities 


led to a small improvement in the forecast. 
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APPENDIX A 


FINITE DIFFERENCE EQUATIONS 


The finite difference equations given here are basically the same 
as those given in Mihok (1974). Additions have been made to include 
special handling of the equations near the pole and the discussion of 
the determination of the gradient of terrain pressure has been changed. 
Friction, diffusion, and heating and moisture source terms are not 
included. In the following equations At is the time increment, AA 
is a five degree longitudinal grid increment and A@ is a 5 degree 
latitudinal grid increment, and Ao = 0.2 is the vertical increment 
between layers. Refer to Figure 1 for the staggered spherical grid 
arrangement and subscripting. Figure 2 illustrates the vertical grid. 
It should be noted that in the following difference equations the k 
subscript on w indicates even sigma levels (0.8, 0.6, 0.4, 0.2) and 
for all other quantities indicates odd sigma levels (0.9, 0.7, 0.5, 
Oeste Oek)\s 


Longitudinal momentum equation 
1. Left hand side 


sue age Oe (1) 
aa 2At 
where 
a 
Re es iat ie 1 tt ee at) 


In the above Equation (1.1) the subscript 
L=i (i even) 


L = i-1 (j odd) 
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2. Right hand side 


a. Term one 


Sa ee 
a cos 6 


oe 4 duvtt cos 8),_ 
= 90 i= 


—_— — 


ull 


Oe) 3 Phas ow Pivtea eto. tae © MH 4 1-1/2, 40k 
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(2) 
where T is computed as in (1.1) and 
— iy 
G arif2,3,k glum) 5 OM a4 1,567, j+1 PT @M 1, 5-1, (2.2) 
O19 4 & is obtained from (2.1) by incremental subscripting. 
vi = (i), +7) +i) +i) are 
i,jtl,k 4 ijk i,j+2,k te ic M, j+1,k : 
vl. . is obtained from (2.2) by incremental subscripting. 
i,j-l,k 
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When computing the above expressions around the 85° and 80° latitude 
circles any term requiring an undefined polar quantity is ignored. 


For example 


fai 7 Cae = 
2 441/2,4,k 


— 


; A . ° * T 
is not included in the computation at 85°S since the On 41/24, 


contains the undefined quantity (um), soy 
> > 


b. Term two 





30 ij 2Ao0 


where the vertical velocity is defined as 


w eget 
ijk Lj sk-L oT ot’ ij 
1 orm 44 4k ak» See ea pea 
a cos®@ 2A 2A6 


and the local chance in terrain pressure as 
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The 7 in the above equations is computed as in (1.1). The w terms 


are computed in the same manner as the 7 terms, 
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When computing w and the local change in terrain pressure at the 
north pole for the thermodynamic and moisture equations the following 


expressions are used, as discussed in Section II-D-2, 





_ ano. 
mole, “polesk=1 re Ge pole * aA036 yw i,35,k! (6) 
and 
on : 2Ac 28 
T = 
GP = - Da ah036 DOM, 35 & (7) 
= =i | 


In order to compute these quantities at the south pole change the sign 
of (v™) and the subscript j from 35 tol. 


c. Term three 


Tuv tan9 (uTTv) 5 4 tand, N 
Ene ee ee (8) 


a a 


3! 


with as in (1.1). 


d. Term four 
Tvf = f (TY) 5 5k (9) 


with T as in (1.1). 


e. Term five 


SS = 
i a cos [RT a at or! = 
ES AE EE 
1 = om, * - ijk i-1,j,k 
A eon Oo Oe ty AN J (10) 


where 7 is obtained from (1.1) and T is computed in the same manner. 


The, 
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Based on Equation (10) of Section II-A the longitudinal derivative 


of terrain height may be expressed 


om, * T 
Gj = [Ad -— Ad_] . (11) 
oA RT2A.—P 0 





In order to obtain ee - Abe] we first use the hypsometric equation 


to form 


- Ad, =e ORT [en T= Qn 1,4! (12) 
The right hand side of (12) can be modified to 
- A$, = RT[2nt,-LnT, + nm gat; 4] (13) 


without changing the equality. 
Bringing T inside the parenthesis and computing it in a special way 


on either side of the velocity point yields 


89 is MO Vat 


Pp ie} 
BU arr kag ay Trront 2" 5 hae (14) 
TL EFT and TRIGHT are special mean temperatures derived over the five 


degree grid distance from the computational velocity point to the 
adjacent mass point either on the right or on the left. The following 
method is used to compute these temperatures, The subscript S indi- 
cates side (either right or left). 

SO See & 2 (nm, ; - £nm.) <r or k=1, (ise. oO = 0.9) then 


meer) as = hey). Cam, - nm.) na 
Ske Sk 2 (ino, )- no.) 
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Equation (15) takes care of the situation where the pressure on the 
sigma surface at the computational point is less than the sigma surface 
pressure at an adjacent mass point. The mean change of temperature 
between the layers is multiplied by the logarithmic slope of the sigma 


surface au - This is used to modify the temperature of the lower 





layer. The end result is to interpolate the value of the geopotential 
at the mass point to the same pressure surface as the computational 
point before integration takes place. 


(2) 235 (nm, , - £nT,) > 0 or k= 5, (0 = 0.1), then 


(Qn, .-2nt..)+ (eno, ~k£n0,_ 4) 


&n(o, - {n o 


_ Cli eee) 
pees? oy So eee 


) ] (16) 
k-1 
Equation (16) comes into play when the pressure on the sigma surface 
at the computational point is less than the pressure on the sigma sur- 
face at an adjacent mass point. The mean change of inberabaee between 
the layers is multiplied by the logarithmic slope of the pressure sur- 


face ape - This is used to modify the upper layer. The end result 





is to interpolate the value of the geopotential at the mass point to the 
same pressure surface as the computational point before integration takes 
place. ' 

This computational scheme for a is intended to reduce truncation 
error when differencing over steep terrain. Even though a flat earth 


is used in this thesis the terrain pressure spatial derivatives using 


this scheme were allowed to remain. 
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Latitudinal momentum equation 
1. Left hand side 


- ntl -,n-l 
out . (Waa ~ (yay 
ae Cae hs 


where 1 is computed as in (1.1). 
2. Right hand side 


a. Term one 


* 1 po civ) _ d(tvv cos 8) 


a cosd, "Bh eee sl 
ees vase a1, jk /2, jk Vague in1 1 M-1/2, 5k 
a eant 2 ZN 
T Eat x eat 
te ee tee ak a, fae ad kote 
‘i 2 2N0 ! 


(18) 
where the un and vm terms are computed as in (2.1) and (2.2). 


b. Term two 


Moa (VIR ty So) ow, . (v..,+.. ) 
d(wv) _ - Wijk Pepe ak “jal jsk-L" ijk 1,j,k-1 
ae aloes 7, ! 2Ac 2Ac ] (19) 


- oT 4 
The w and GP terms are computed in the same manner as they are 
for (3), the vertical acceleration term in the latitudinal momentum 
equation. 


c. Term three 


= (umu).., tan 6. 
zs ae Xe —* j (20) 


where 1 is computed as in (1.1). 
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d. Term four 


juf = £5 (Tu) 5 5, (21) 
where 7 is obtained as in (1.1). 
e. Term five 
i OTT O06, - 
— ql RT 36 + 17 Pye = 
ae ae ) 
do ame | = L,jtl,k _“L,j-1,k 
SOE a Gel et My IN) ] (22) 
where Lei (j even) 


L=i-l (j odd) 


T is obtained from (1.1) and T is computed in a similar manner, 
The computation of eu ny is carried out in the same manner as the com- 
putation of a * for Equation (10) with the exception that the Cl 
in this case is computed along a meridian instead of a latitude circle. 
Refer to Equations (10), (11), (12), (13), (14), (15), and (16). 
Thermodynamic equation 

1. Left hand side 


pe n-1 

ee) oe Leite a 14k es 

ae Tae 
This term is computed directly at the poles. 

2. Right hand side 

a. Term one 

Set ed acco ey Os ie 
a cos @ or Pls) = S aldés es 2Ar 
2 (Tvt 7), ‘41 4.0088 54 = ren), 0089; 1, As 
a 


79 





mae sadsate e nt betrrtore’ 

ates Gilt. ct ee Omh' aft nut 

an | : 

€ ONS Sad wort pce’) os 
hate Maa 


~Oe0e% obera lint a. Hagan! 


tts) 


; i Ma oe 
i 
’ dunit ; 
i 4 
eat) ; 
La ie 
se 
¥ ea 
ij : 
, 
; 4 
{ reves 
‘. 
) 
ty ’ 
7ou3 
; 
? 
me 





where 
L=i (j odd) 


L = itl (j even) 


T is computed as in (1.1) and T is computed in a similar manner. At 


85°S (vm), ghthyk is undefined and not included in the difference equa- 
> > 


tion. At 85°N (Tvm) is undefined and not included. The 


L, j+1,k 
computation of (24) at the poles is carried out as discussed in Section 


II-D-2 (Polar Finite Differencing). At the north pole the expression 


used for term one, (24), is 


36 
2 2 == 
R.H.S. term 1 = - Saye be (IVT) 5 35 & (25) 
i=1 
In order to compute this term at the south pole change the sign of 
(Ivt) and the subscript j to 1. 


b. Term two 
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where Equations (4) and (5) are used to compute w. This term is 
computed directly at the poles. 


c. Term three 


oT 1 oT oT = 
rac] ea + ry eye, + v Gp) cos 01)] = 


* Tcosd 2X 2X 
hee COT) 7 jet VOL j-1,k _ nae M1, f+ YL, j-1, 0, 
2A0 ij 2h0 
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where 
b= i (odd) 


L = itl (j even). 


w is computed by (4), a by (5), and tw by (1.1). Note a special 


scheme is used to compute the derivatives of terrain pressure in this 


equation in order to reduce truncation error. At 85°S - vm quik and 
> > 
“vy ee. are undefined and not included in the difference equation. 
> > 
° — ; ° 
At 85°N Wily etd and Bat elr are undefined and not included. The 


computation of (27) at the poles is carried out as discussed in Section 
II-D-2 (Polar Finite Differencing). At the north pole the expression 


used for term three (27) is 


R.H.S. term 3 = Rtpole,k |-" Sled | “pole kl "pole 





CO, 2 
36 
oT 2 a 
* a, (#2 + 5A036 De LOE) saa > r,35%1,35,4)| £28) 
pole 
i=1 | 


In order to compute this term at the south pole change the sign of 
(vt) and (mv) and change the subscript j to l. 


Moisture equation 


1. Left hand side 


n+1 n-1 
dqm _ (47) 5 5% i CULE (29) 
oe: 2 2At 


“This term is computed directly at the poles. 
2. Right hand side 


a. Term one 
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(2ounD). 4 d(vTg cos Bi stat 1 
a cos 6 00 


L, j+1,k 
+ DAO (30) 


(qvt) cos®. |) - (qvT) | pa ORO. 
where 
y= et (j odd) 


L = i+l (j even) 


tT is computed as in (1.1) and q is computed in a similar manner. 


At 85°S (qv) is undefined and at 85°N (qvm) is 


L,j-1,k L,jtl,k 


undefined. Undefined terms are not included in the difference equations. 
The computation of (30) at the poles is as discussed in Section II-D-2 
(Polar Finite Differencing). At the north pole the expression used for 
term one, (30), is 

36 


es 2 » aC 


i=1 
In order to compute this term at the south pole reverse the sign of 


(qv) and change the subscript j to l. 


b. Term two 





oad), rf aie poke *Si ji? 7 Wing ker Gaje’ Gi, jue 
30 


ij 2Ao0 ] (32) 


where Equations (4) and (5) are used to compute w. This term is com- 
puted directly at the poles. 
Hydrostatic equation 

The hydrostatic equation is integrated to form the hypsometric 


equation te 
Ad = -RT A &n 0 (33) 
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This is converted to a difference equation which is a modified form of 


that used by Kesel and Winninghoff (1972). 


4 4,ke1 ” 


T, 54, (2n9 +knm, .)+T, (2no, tent, .) 


0. 
k+]l, i k j,k+1 
eggs AT Coane) Uf ] (34) 
ijk O £no). + kno 47 “F zis, 
For the lowest layer (34) is modified to form 
os51 = SBS &no, (35) 


Both (34) and (35) are computed directly at the poles. 
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APPENDIX B 


TERMS USING FOURTH ORDER FINITE DIFFERENCING 
IN MIXED SECOND AND FOURTH ORDER SCHEME 


Williams (1972) has concluded that fourth order differencing used 
only in the advection terms and the divergence term of the continuity 
equation will improve the phase speed of meteorological waves while 
keeping the necessary reduction in the time step to a minimum. The 
following equations are the advection terms and the divergence term in 
fourth order finite difference form as developed by Mihok (1974) for 
use in the global model. 


The advection term for the latitudinal momentum equation is 


P 1 yarn heat 1 "G41 /2.4,k dak t-1,4,k0  "i-1/2.4.k 
a cos 8, 3 2 * 2Ar 


2 * 4A 


(u;., tu; 44 OT; +49 geOs9 a o> CU hy +), WY" 9% 
+. 2+ 4A0 


(1) 


where the ut ig vt and vt, are 


tJ 24k eke, 4 i,jtl,k foil, k 


given by (2.1) and (2.2) in Appendix A. Note that 
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ul ‘ is obtained from (1.1) by incremental subscripting. 
i-1,j,k 


—_ 


atk 


ae a ay = 
Be EE Og 141, 542,62 544, u! 2) 


and vi; 42k is obtained from (1.2) by incremental subscripting. 
> > 
If one substitutes v everywhere for u in (1), the result is the 


longitudinal momentum equation. 


The advection term for the thermodynamic equation is given by 


1 4 (Tum) 544 pk 7 TUM ay 

of ee Seis Nl EA Cee ee 13S 

a cos = 5 { 2Ar 

Z (Ivt) ; a qo089 547 - (Tum) . +4 rae et 
2A8 


—-*k —_% —-*k —_* 
ee ele atte AOE) a 
3 


ak osek —*) —_% 
(T [vt] ); .,5 ,cos®.,, - @ [vt] ), ._, ,cos0._, 
ee a cL (2) 


‘where 7 is computed as in (1.1) of Appendix A and T is computed in 


a similar manner. Note that 
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7 liz +T +7 +T ] (757) 
i-1,j,k 4 N,j,k N,jtl,k 0,j,k N,j,k y 
7" aig +T + 7 +7 
Miso. WET Geek LS OM, F2,k * Spe g41, ic! (2.3) 
and T is obtained from (2.3) by incremental subscripting. 


i,j-2,k 


In the above equations 
L = itl (j odd) 
L = it+2 (j even) 
M = i+2 (j odd) 
M = i+2 (j even) 
N = i-l (j odd) 
N.= i (j even) 
O=i (j odd) 
O = i-l (j even) 
The expressions for [ur] can be obtained by substituting uT in 
for T in the above equations. : 
Quantities with a bar, e.g. T, are applicable at velocity points 
while starred quantities, e.g. T » are applicable at mass points. The 


advection term for the moisture equation is the same as (2) with q 


substituted for T. 
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The divergence term of the continuity equation is 


— 


1 4 OM in dk -@T), fs vt), 41 p09 54 ro), 4cos8. -1 
a eo 3 2Ar aS 
—_* oa ae 8 ag ‘ 
Die trig aed  tetatale oe Wl Mica OP ens tea) go oO 9 
3 4A 4A6 
e (3) 


where 7 is given by (1.1) in Appendix A and the [ur]* quantities 

are computed as in the thermodynamic equation. } 
Mihok (1974) blended the mixed second and fourth order scheme with 

the normal second order scheme at 75° latitude. This was done to avoid 


differencing over the poles. 
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